Data Analysis Notebook

#Load in MBB - Statistic Tracking Report csv
stats <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/MBB - Statistic Tracking Report.csv")

# Load in ACWR - Kinexon MBB
kinexon_ACWR <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/ACWR - Kinexon MBB.csv")

#Loading in VALD - Dynamo MBB
VALD_dynamo <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/VALD - Dynamo MBB.csv")

# Load in VALD - ForceDecks MBB csv
VALD_ForceDecks <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/VALD - ForceDecks MBB.csv")

# Load in Kinexon Session
kinexon_session <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/Kinexon Session MBB.csv")

Data Cleaning

# Remove columns where all data points are NA
stats_clean <- stats |>
  select(where(~ !all(is.na(.))))


# If all three date columns are equal, remove the redundant ones
if (all(stats_clean$Date == stats_clean$Date.Reverse, na.rm = TRUE) &&
    all(stats_clean$Date == stats_clean$Session.Date, na.rm = TRUE)) {
  
  stats_clean <- stats_clean |>
    select(-c(Date.Reverse, Session.Date))
}


# Convert if Date column is NOT already in Date format
if (!inherits(stats_clean$Date, "Date")) {
  stats_clean <- stats_clean %>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(stats_clean$Date)

stats_clean <- stats_clean[!is.na(stats_clean$Plus..Minus),]

# Add more stats
stats_clean <- stats_clean %>%
  filter(Date >= as.Date("2024-09-20", '%Y-%m-%d')) %>%
  group_by(anon_id) %>%
  mutate(Player.Average.Rebounds = mean(na.omit(Rebounds)),
         Player.Average.Steals = mean(na.omit(Steals)),
         Player.Average.Assists = mean(na.omit(Assists)),
         Player.Average.Turnovers = mean(na.omit(Turnovers)),
         Player.Average.Points.Scored = mean(na.omit(Points.Scored)),
         Player.Median.Rebounds = median(na.omit(Rebounds)),
         Player.Median.Steals = median(na.omit(Steals)),
         Player.Median.Assists = median(na.omit(Assists)),
         Player.Median.Turnovers = median(na.omit(Turnovers)),
         Player.Median.Points.Scored = median(na.omit(Points.Scored))) %>%
  ungroup() %>%
  group_by(Date) %>%
  mutate(Team.Rebounds = sum(na.omit(Rebounds)),
         Team.Steals = sum(na.omit(Steals)),
         Team.Assists = sum(na.omit(Assists)),
         Team.Turnovers = sum(na.omit(Turnovers)),
         Team.Points.Scored = as.numeric(str_split(Game.Score, "-")[[1]][1])) %>%
  ungroup() %>%
  mutate(Team.Average.Rebounds = mean(Team.Rebounds),
         Team.Average.Steals = mean(Team.Steals),
         Team.Average.Assists = mean(Team.Assists),
         Team.Average.Turnovers = mean(Team.Turnovers),
         Team.Average.Points.Scored = mean(Team.Points.Scored)) %>%
  mutate(Player.Good.Game = ifelse(Rebounds >= Player.Median.Rebounds & Steals >= Player.Median.Steals & Assists >= Player.Median.Assists,1,0),
         Player.Bad.Game = ifelse(Rebounds < Player.Median.Rebounds & Steals < Player.Median.Steals & Assists < Player.Median.Assists,1,0),
         Team.Good.Game = ifelse(Team.Rebounds >= median(Team.Rebounds) & Team.Steals >= median(Team.Steals) & Team.Assists >= median(Team.Assists),1,0),
         Team.Bad.Game = ifelse(Team.Rebounds < median(Team.Rebounds) & Team.Steals < median(Team.Steals) & Team.Assists < median(Team.Assists),1,0)) %>%
  select(-c(Sport, Day.of.Week, Month, Month.with.Number, Year, Number.of.Competitions.Won, Total.Daily.Competitions.Available, Total.Competitions.Won, Total.Game.Minutes, Offensive.Rebounds, Defensive.Rebounds, Buff.Plays, Percentage.Daily.of.Competitions.Won)) %>%
  filter(!(anon_id %in% c("ID_10", "ID_14", "ID_15", "ID_18", "ID_54", "ID_64")))

# Remove columns where all data points are NA
kinexon_ACWR_clean <- kinexon_ACWR |>
  select(where(~ !all(is.na(.))))



# If all three date columns are equal, remove the redundant ones
if (all(kinexon_ACWR_clean$Date == kinexon_ACWR_clean$Date.Reverse, na.rm = TRUE) &&
    all(kinexon_ACWR_clean$Date == kinexon_ACWR_clean$Session.Date, na.rm = TRUE)) {
  
  kinexon_ACWR_clean <- kinexon_ACWR_clean |>
    select(-c(Date.Reverse, Session.Date))
}


# Convert if Date column is NOT already in Date format
if (!inherits(kinexon_ACWR_clean$Date, "Date")) {
  kinexon_ACWR_clean <- kinexon_ACWR_clean %>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(kinexon_ACWR_clean$Date)


# Selecting only relevant columns
kinexon_ACWR_clean <- kinexon_ACWR_clean %>%
  select(anon_id, Date, Day.of.the.Week, Month, Year, Position, Session.Type, Session.Description, Practice.Classification, Last.7.Day.Acceleration.Load, Today.s.Acceleration.Load, Acceleration.Load.ACWR, Acceleration.Load.EWMA.ACWR, Acute.Total.Acceleration.Load, Acute.Total.Acceleration.Load.EWMA, Chronic.Total.Acceleration.Load, Chronic.Total.Acceleration.Load.EWMA, Total.Low.Acceleration.Load, Total.Medium.Acceleration.Load, Total.High.Acceleration.Load, Total.Very.High.Acceleration.Load, Last.7.Day.Total.Jump.Load, Today.s.Jump.Load, Total.Jump.Load.ACWR, Acute.Total.Jump.Load, Acute.Total.Jump.Load.EWMA, Chronic.Total.Jump.Load, Chronic.Total.Jump.Load.EWMA, Last.7.Day.Total.Minutes, Today.s.Minutes, Minutes.ACWR, Minutes.EWMA.ACWR, Acute.Minutes, Acute.Minutes.EWMA, Chronic.Minutes, Chronic.Minutes.EWMA, Total.Distance.EWMA.ACWR)

#selecting important variables, removing rows with missing data, only selecting hand movements
VALD_dynamo_clean <- VALD_dynamo %>%
  filter(Body.Region == "Hand") %>%
  mutate(Right.Side = ifelse(Repetition.Type.Laterality=="RightSide",1,0)) %>%
  select(anon_id,   Date, Max.Force.Newtons, Avg.Force.Newtons, Max.Impulse.Newton.Seconds,  Max.Rate.Of.Force.Development.Newtons.Per.Second,  Avg.Rate.Of.Force.Development.Newtons.Per.Second,   Min.Time.To.Peak.Force.Seconds, Start.Offset.Seconds, Repetition.Duration.Seconds,  Repetition.Max.Force.Newtons,   Impulse.Newton.Seconds, Rate.Of.Force.Development.Newtons.Per.Second,   Time.To.Peak.Force.Seconds,Right.Side) %>%
  na.omit()

# Convert if Date column is NOT already in Date format
if (!inherits(VALD_dynamo_clean$Date, "Date")) {
  VALD_dynamo_clean<- VALD_dynamo_clean%>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(VALD_dynamo_clean$Date)


#This data set contains information from each individual hand in a strength test meaning that there are two measurements per day for each player. This makes it hard to understand in plots especially when plotted like time series data. For future reference, it might be useful to take an average of both sides or note that there is a difference in each hand and plot the left and rights sides separately.

#### Note that you have to group on id and the date in order to get a valid average of both sides for a specific day for each individual player. 

VALD_ForceDecks_clean <- VALD_ForceDecks %>%
  select(where(~ !all(is.na(.)))) %>%
  mutate(X=1:3384) %>%
  filter(Test.Type == "DJ") %>%
  select(X,anon_id, Date, Week, Position, Test.Type, Weight..kg., Trial, Flight.Time..s., Jump.Height..Flight.Time...cm., Eccentric.Concentric.Mean.Force.Ratio, Contact.Time..s., RSI..JH..Flight.Time..Contact.Time., RSI..Flight.Time.Contact.Time., Drop.Height..cm.)%>%
  mutate(RSI..m.per.s. = (Jump.Height..Flight.Time...cm./100)/Contact.Time..s.,
         Date = as.Date(Date, '%m/%d/%Y')) %>%
  group_by(Date, anon_id) %>%
  mutate(Weight..kg. = mean(na.omit(Weight..kg.)),
         Player.Average.RSI..m.per.s. = mean(RSI..m.per.s.)) %>%
  ungroup() %>%
  group_by(Date) %>%
  mutate(Team.Average.RSI..m.per.s. = mean(na.omit(RSI..m.per.s.)))

# Remove columns where all data points are NA
kinexon_session_clean <- kinexon_session |>
  select(where(~ !all(is.na(.))))


# Convert if Date column is NOT already in Date format
if (!inherits(kinexon_session_clean$Date, "Date")) {
  kinexon_session_clean <- kinexon_session_clean %>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(kinexon_session_clean$Date)

[1] "Date"

# Filter to only include data from this season
kinexon_session_clean <- kinexon_session_clean |>
  filter(Date >= as.Date("2024-6-01"))

##Analysis

Question 1

Are changes in RSI related to team game performance?

Team Wins and Losses

#looking at the average RSI as well as the adjusted variance
#functions to calculate mean and variance of time series found in "Time-Series-Functions.Rmd"
team_mean <- mean(current_season$Team.Average.RSI..m.per.s.)

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'mean': object 'current_season' not found

The plot above shows the team’s average RSI per week throughout the 2024-2025 season. This plot shows that there are large variations in team’s RSI values throughout the season but there does not seem to be any statistically significant underlying trend in the data. The lines in the plot represent the mean and one standard deviation from the mean throughout the season. As you can see, there is a lot of variability throughout the data. The patterns that appear though could potentially suggest an underlying ARIMA or other model for the data.

#calculating mean RSI values for each day and attaching them to correct dates
Team.means <- unique(current_season$Team.Average.RSI..m.per.s.)
Dates <- unique(current_season$Date)
IDs <- unique(current_season$anon_id)
Team.Change <- rep(0, 25)
current_season$Change.Team.Average.RSI..m.per.s. <- rep(NA, 580)

#calculating the difference in RSI from the previous measurement and attaching as a new column in the current season data set
for(i in 2:25){
  Team.Change[i] = Team.means[i] - Team.means[i+1]
  current_season[current_season$Date == Dates[i],19]<- Team.Change[i]
}

#making binary column of whether RSI increased or not
current_season <- current_season %>%
  mutate(Team.RSI.Increase = ifelse(Change.Team.Average.RSI..m.per.s. > 0, 1, 0))


#calculating the difference in RSI from the previous measurement for each player and attaching as a new column in the current season data set
current_season$Change.Player.Average.RSI..m.per.s. <- rep(0,580) #making new column

for(i in 1:11){ #cycle through players, calculate difference from last measurement
  Dates <- unique(current_season[current_season$anon_id == IDs[i],]$Date)
  Player.Means <- unique(current_season[current_season$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.) #taking measurements throughout season
  Player.Change <- rep(0,length(Player.Means))
  ifelse(i != 14, Player.Change <- rep(0,length(Player.Means)), 0) #don't need this line any more but I am not going to delete it

  for(j in 1:length(Player.Change)){ #calculates changes from last measurements
    Player.Change[j] <- Player.Means[j] - Player.Means[j+1]
    }
  #attaching RSI changes in new column  
  for(k in 1:length(Player.Change)){
    current_season[current_season$Date == Dates[k] & current_season$anon_id == IDs[i],]$Change.Player.Average.RSI..m.per.s. = Player.Change[k]
  }
}

#making indicator variable for whether RSI increased from last measurement or not
current_season <- current_season %>%
  mutate(Player.RSI.Increase = ifelse(Change.Player.Average.RSI..m.per.s. > 0, 1, 0))

Maybe Delete

#represents games throughout the season, will merge into one data set
Win <- c(1,1,1,1,1,0,1,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,1,1,0,0)
All_Dates <- as.Date(c("2024-10-19","2024-11-4", "2024-11-8", "2024-11-13", "2024-11-17", "2024-11-25","2024-11-26", "2024-11-27", "2024-12-2", "2024-12-7", "2024-12-13", "2024-12-21", "2024-12-30", "2025-1-4", "2025-1-8", "2025-1-12", "2025-1-15", "2025-1-18", "2025-1-21", "2025-1-25", "2025-1-28", "2025-2-2", "2025-2-5", "2025-2-8", "2025-2-11", "2025-2-15", "2025-2-18", "2025-2-22", "2025-2-24", "2025-3-2", "2025-3-5", "2025-3-8", "2025-3-11","2025-3-12", "2025-3-13", "2025-4-1"),"%Y-%m-%d")
Regular <- as.Date(c("2024-10-19","2024-11-4", "2024-11-8", "2024-11-13", "2024-11-17", "2024-12-2", "2024-12-7", "2024-12-13", "2024-12-21", "2024-12-30", "2025-1-4", "2025-1-8", "2025-1-12", "2025-1-15", "2025-1-18", "2025-1-21", "2025-1-25", "2025-1-28", "2025-2-2", "2025-2-5", "2025-2-8", "2025-2-11", "2025-2-15", "2025-2-18", "2025-2-22", "2025-2-24", "2025-3-2", "2025-3-5", "2025-3-8"),"%Y-%m-%d")
Tournament <- as.Date(c("2024-11-25","2024-11-26", "2024-11-27"),"%Y-%m-%d")
Post_Season <- as.Date(c("2025-3-11","2025-3-12", "2025-3-13", "2025-4-1"),"%Y-%m-%d")

#summarizes all games throughout the season
record <- data.frame(Date = All_Dates, Win = Win) %>%
  mutate(Type = case_when(
    Date %in% Regular ~ "Regular",
    Date %in% Tournament ~ "Tournament",
    Date %in% Post_Season ~ "Post-Season",
    TRUE ~ "NA"
  ))

#plotting regular game outcomes against RSI values
ggplot(current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line(linewidth = 0.75) +
  geom_vline(xintercept = record$Date[record$Type == "Regular"], linewidth = 0.2, color = ifelse(record$Win[record$Type == "Regular"] == 1, "palegreen3", "red3"))

#plotting tournaments or post season games with average team RSI values
ggplot(current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line(linewidth = 0.75) +
  geom_vline(xintercept = record$Date[record$Type %in% c("Tournament","Post-Season")], linewidth = 0.2, color = ifelse(record$Win[record$Type %in% c("Tournament", "Post-Season")] == 1, "palegreen3", "red3"))

These plots show that there may not be much predictive power for team game performance from team RSI values themselves. Instead, it may be useful to look at the changes in RSI before each game or look at a normalized plot of the team’s average RSI to get a look at what it would look like without as much noise.

The data was normalized with a moving average process which removed some of the noise present in the data. For this plot, the green lines represent a game in which the team’s stats were higher than the median for steals, assists, and rebounds. Red lines represent games where the total team sum of steals, assists, and rebounds are lower than the season median. Again, There does not seem to be a logical association between team average RSI and in game statistics. Instead, we find that games that had better in game performance are associated with lower RSI recordings before the game and games that had worse in game performance are associated with higher RSI recording just before the game. It’s important to note that there gaps between measurements that could have missed important spikes or dips in RSI for the team that would give us a more clear understanding of the relationship. Previous research that was able to identify significant relationships between in game performance and RSI measurements were able to collect it before and after each event.

#plotting the change in team average RSI throughout the season
ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)


#plotting the changes in team average RSI throughout the season with green games being games where the team performed better than their median performance in explosiveness metrics and red being games where the team performed below median for explosiveness metrics
ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))
#  geom_text(data = stats_clean, x = Date, y = 0.19, label = Date, angle = 90,size = 2.75)

This plot shows the change in team average RSI measurements throughout the season. Each data point represents how much the average RSI changed from the previous measurement. This plot helps to make sense of how the RSI of the team changed throughout the season. We can see that while increases were more frequent, they were smaller than the drops in RSI that occurred throughout the season. With this though, we can see that there is potentially a relationship between changes in RSI and in team game performance. This is because visually, it seems that games where the team had a “good” game tend to have an increase in RSI measured right before the game and games that are “bad” tend to have a decrease in RSI right before the game or a less dramatic increase in RSI than games that are “good”.

Maybe Delete

ggplot(data = current_season, aes(x=Date, y=Team.Average.RSI..m.per.s.)) +
  geom_line(aes(color = Team.RSI.Increase)) +
  geom_vline(xintercept = record$Date[record$Win == 1], linewidth = 0.2, color = "palegreen3")

ggplot(data = current_season, aes(x=Date, y=Team.Average.RSI..m.per.s.)) +
  geom_line(aes(color = Team.RSI.Increase), linewidth = 0.75) +
  geom_vline(xintercept = record$Date[record$Win == 0], linewidth = 0.2, color = "red3")

Team Explosiveness Metrics

ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Points.Scored > median(stats_clean$Team.Points.Scored), "palegreen3", "red3"), linewidth = 0.2)


ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0) +
   geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Points.Scored > median(stats_clean$Team.Points.Scored), "palegreen3", "red3"), linewidth = 0.2)

This plot represents if the team scored higher or lower than the median team score in that game. Green represents higher and red represents lower than the median. This plot suggests that lower than median scores are associated with drops in RSI before the game (found in the middle of the season). As the season goes on, we can see that even though the team’s RSI goes up, the team has a consistently lower than median score. This could be indicative of trouble recovering from fatigue found in earlier research.

#Looking at games that had team total rebounds above or below median team season value
ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Rebounds > median(stats_clean$Team.Rebounds), "palegreen3", "red3"), linewidth = 0.2)

ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)  +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Rebounds > median(stats_clean$Team.Rebounds), "palegreen3", "red3"), linewidth = 0.2)

#Looking at games that had team total steals above or below median team season value
ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Steals > median(stats_clean$Team.Steals), "palegreen3", "red3"), linewidth = 0.2)

  ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)  +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Steals > median(stats_clean$Team.Steals), "palegreen3", "red3"), linewidth = 0.2)

#Looking at games that had team total assists above or below median team season value
ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Assists > median(stats_clean$Team.Assists), "palegreen3", "red3"), linewidth = 0.2)

ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)  +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Assists > median(stats_clean$Team.Assists), "palegreen3", "red3"), linewidth = 0.2)

ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))

Applying regularization techniques to team data

Dates <- unique(current_season$Date)
#smooth <- exponential_smoothing(unique(current_season$Team.Average.RSI..m.per.s.),2,0)
ma <- moving_average(unique(current_season$Team.Average.RSI..m.per.s.), 1)
#x53 <- five_three_x(unique(current_season$Team.Average.RSI..m.per.s.))

#current_season$Team.Average.RSI.Smooth <- rep(NA,580)
current_season$Team.Average.RSI.Moving.Average <- rep(NA, 580)
#current_season$Team.Average.RSI.53X <- rep(NA, 580)

for(i in 1:25){
  #current_season[current_season$Date == Dates[i],21]<- smooth[i]
  current_season[current_season$Date == Dates[i],23]<- ma[i]
  #current_season[current_season$Date == Dates[i],23]<- x53[i]
}

ggplot(data=current_season, aes(Date, Team.Average.RSI.Moving.Average)) +
  geom_line() +
geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))

Assessing effectiveness of RSI as predictor of in game success

#getting dates that had games on them
Game.Days <- unique(stats_clean$Date)

#numbering which game in the season it was
Game.Number <- seq(30, 1, -1)

#making empty columns to be filled with game number
stats_clean$Game.Number <- rep(NA, 315)
current_season$Game.Number <- rep(NA, 580)

#going through games data set and labeling which game was which throughout the season
## first, second, third, ... 30th
for(i in 1:30){
  stats_clean[stats_clean$Date == Game.Days[i],]$Game.Number <- Game.Number[i]
}

#going through practices data set and labeling the closest game that the practice was before
## will be used for grouping games together later
for(i in 1:30){
  current_season[current_season$Date < Game.Days[i],]$Game.Number = 31-i
}

#making one data set with everything together to assess effectiveness of RSI as a predictor
## grouping by the game played or led up to and the player id
Q1_data <- left_join(x=current_season, y=stats_clean, by=c("anon_id", "Game.Number"))
Q1_data <- Q1_data %>%
  filter(Team.Good.Game %in% c(0,1))
Q1_data <- Q1_data[13:540,]
Q1_data <- Q1_data[-c(468,469,514,515,516),]

#set seed to make reproducible, splitting data into training and testing sets (75-25 split)
set.seed(123)
rows <- sample(1:nrow(Q1_data), nrow(Q1_data)*0.75, replace=FALSE)
Q1_train <- Q1_data[rows,]
Q1_test <- Q1_data[-rows,]

#model to predict whether game was a good game or not based on team average RSI increase leading up to that game
Q1_model_good_game <- glm(Team.Good.Game~Change.Team.Average.RSI..m.per.s., data=Q1_train, family="binomial")
summary(Q1_model_good_game)

#using the model above to make prediction on whether it was a good game or not and calculating CER (using 0.5 as threshold, may adjust later)
### Changed threshold to 0.43 based off of ROC curve done below
Q1_test$Team.Good.Pred <- ifelse(predict.glm(Q1_model_good_game, Q1_test, type="response") >0.46,1,0)

Q1_test %>%
  summarize(CER = mean(na.omit(Team.Good.Pred != Team.Good.Game)))

#building binomial model to predict whether it was a bad team game or not
Q1_model_bad_game <- glm(Team.Bad.Game~Change.Team.Average.RSI..m.per.s., data = Q1_train, family = "binomial")
summary(Q1_model_bad_game)

#making prediction based on the model and calculating the CER for the model on test data
Q1_test$Team.Bad.Pred <- ifelse(predict.glm(Q1_model_bad_game, Q1_test, type="response") >0.05,1,0)

Q1_test %>%
  summarize(CER = mean(na.omit(Team.Bad.Pred != Team.Bad.Game)))

Cross Validating the Team Models and Calculating Sensitivity and Specificity Rates w/ ROC

set.seed(123)
#model to be used for cross validation and ROC curve
model_good_team_cv <- glm(Team.Good.Game ~ Change.Team.Average.RSI..m.per.s., data=Q1_data, family="binomial")

#Making predictions for ROC curve
#$Team.Good.Game.Pred <- predict(model_good_team_cv, type="response")

#making ROC curve
#pred <- prediction(predictions = Q1_data$Team.Good.Game.Pred,
#                   labels = Q1_data$Team.Good.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)


#cost function for CER in cross validation
cost <- function(obs, pred){
  mean((pred <= 0.46) & obs==1 | (pred > 0.46) & obs==0)
}

#cross validating model
good_team_cv <- cv.glm(data=Q1_data,glmfit=model_good_team_cv,cost,K=5)

#extract average error
good_team_cv$delta[1]

Q1_data$Team.Good.Game.Pred <- ifelse(predict(model_good_team_cv, type="response")>0.46,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Team.Good.Game.Pred, y=Q1_data$Team.Good.Game)

set.seed(123)
#model to be used for cross validation and ROC curve
model_bad_team_cv <- glm(Team.Bad.Game ~ Change.Team.Average.RSI..m.per.s., data=Q1_data, family="binomial")

#Making predictions for ROC curve
#Q1_data$Team.Bad.Game.Pred <- predict(model_bad_team_cv, type="response")

#making ROC curve
#pred <- prediction(predictions = Q1_data$Team.Bad.Game.Pred,
#                   labels = Q1_data$Team.Bad.Game)
#perf <- performance(pred,measure='fnr',x.measure='fpr')
#plot(perf, colorize=TRUE)


#cost function for CER in cross validation
cost <- function(obs, pred){
  mean((pred <= 0.05) & obs==1 | (pred > 0.05) & obs==0)
}

#cross validating model
bad_team_cv <- cv.glm(data=Q1_data,glmfit=model_bad_team_cv,cost,K=5)

#extract average error
bad_team_cv$delta[1]

Q1_data$Team.Bad.Game.Pred <- ifelse(predict(model_bad_team_cv, type="response")>0.05,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Team.Bad.Game.Pred, y=Q1_data$Team.Bad.Game)

In order to limit variability and mimic past research, binary variables were created. Other research used top running speed in games to measure explosiveness of an athlete at a particular time and were compared to median top in game running speed throughout the season to define success. The data sets given did not have in game running speed so a combination of steals, assists, and rebounds were used as a proxy for a measure of explosiveness. A binary variable was created that measured if the team performed above or below median in all three metrics (1) or if one or more of the listed metrics were below the season team median (0). This is what will be predicted through changes in RSI of the team throughout the season. Shown above are the cross validated predicted error rates for the models that are used to predict whether the team will play above their median game play or below their median game play in a given game.

For the first model, it is trying to predict whether or not the team will have their steals, rebounds, and assists, all above their season median for a given game where the only predictor is the change in RSI measured the previous day. This model has a cross validated classification error rate of around 0.4 meaning that it will make a wrong prediction around 40% of the time. Out of the misclassifications, the model will generally tend to predict that a game will not be above median when it is. In fact, the model will predict that a game will not be above median when it actually is, about 3 times as often as it will correctly predict an above median game will be above median. Due to this major lack of sensitivity, it’s somewhat safe to say that changes in RSI may not be a good predictor for in game success based on the constraints made in this analysis.

The second model is trying to predict whether or not the team will have the total rebounds, assists, and steals for the game all below median or not. With this model, the cross validated classification error rate is around 0.2. This suggests that the model will make an incorrect classification around 20% of the time. Out of the incorrect classifications all of them were times in which in the model predicted for the team to perform below their median in explosiveness metrics in the game but the team performed above their median in at least one metric. This suggests that this model is making predictions that the team will perform “poorly” or below median more frequently than what is actually happening. In fact, when the model predicts that the team will perform below their season median game play, the predictions will be incorrect more often than correct. It’s important to note that this model correctly identified all of the games in which the team performed below median. This suggests that there is predictive power in changes in RSI for in game performance that is measured to be below median for the team.

For the data points that were misclassified by the model, it may be useful to look at the metrics for those games and identify how the team performed. Since this analysis defined a bad game as a game in which the team performed below median in steals, assists, and rebounds or not, it may be useful to see which metrics they performed below median in games that were predicted to be “bad games” but where not below median in all three metrics. The model could have potentially identified games that were poor but were not poor enough per say to be classified as a “bad game” in this analysis.

Question 2

Are changes in RSI related to individual statistical game performance?

#Looking at games that had above or below median rebounds throughout the season
ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42"], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42",]$Rebounds >= median(stats_clean[stats_clean$anon_id == "ID_42",]$Rebounds), "palegreen3", "red3"))

#Looking at games that had above or below median steals throughout the season
ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42"], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42",]$Steals >= median(stats_clean[stats_clean$anon_id == "ID_42",]$Steals), "palegreen3", "red3"))

#Looking at games that had above or below median assists throughout the season
ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42"], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42",]$Assists >= median(stats_clean[stats_clean$anon_id == "ID_42",]$Assists), "palegreen3", "red3"))

ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1)], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1),]$Player.Good.Game == 1, "palegreen3","red3"))

ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Change.Player.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept=0) +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1)], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1),]$Player.Good.Game == 1, "palegreen3","red3"))

Applying regularization techniques to player data

smooth <- exponential_smoothing(unique(current_season$Team.Average.RSI..m.per.s.),2,0)
ma <- moving_average(unique(current_season$Team.Average.RSI..m.per.s.), 1)
x53 <- five_three_x(unique(current_season$Team.Average.RSI..m.per.s.))

current_season$Team.Average.RSI.Smooth <- rep(NA,580)
current_season$Team.Average.RSI.Moving.Average <- rep(NA, 580)
current_season$Team.Average.RSI.53X <- rep(NA, 580)

for(i in 1:25){
  current_season[current_season$Date == Dates[i],21]<- smooth[i]
  current_season[current_season$Date == Dates[i],22]<- ma[i]
  current_season[current_season$Date == Dates[i],23]<- x53[i]
}

ggplot(data=current_season, aes(Date, Team.Average.RSI.Moving.Average)) +
  geom_line() +
geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))

Assessing effectiveness of RSI as predictor of in game success

Good Games Models

#building completely pooled model for player RSI changes and good games
Q2_model_good_game_pooled <- glm(Player.Good.Game~Change.Player.Average.RSI..m.per.s., data=Q1_train, family = "binomial")
summary(Q2_model_good_game_pooled)

#making predictions with the pooled model for good games
Q1_test$Player.Good.Pooled.Pred <- ifelse(predict(Q2_model_good_game_pooled, Q1_test)>0.5,1,0)

#calculating the CER of the model on the test data
#Q1_test %>%
# summarize(CER = mean(Player.Good.Game != Player.Good.Pooled.Pred))

#making unpooled model for players
Q2_model_good_game_unpooled <-glmer(Player.Good.Game~Change.Player.Average.RSI..m.per.s. + (Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_train, family = "binomial")
summary(Q2_model_good_game_unpooled)

#making predictions for unpooled model on the test data (calculated with probability)
Q1_test$Player.Good.Unpooled.Pred <- ifelse(predict(Q2_model_good_game_unpooled, Q1_test, type="response")>0.46,1,0)

#calculating CER of unpooled model
Q1_test %>%
  summarize(CER = mean(Player.Good.Unpooled.Pred != Player.Good.Game))

#building completely pooled model for player RSI changes and bad games
Q2_model_bad_game_pooled <- glm(Player.Bad.Game~Change.Player.Average.RSI..m.per.s., data=Q1_train, family = "binomial")
summary(Q2_model_bad_game_pooled)

#making predictions for bad games with pooled model on test data
Q1_test$Player.Bad.Pooled.Pred <- ifelse(predict(Q2_model_bad_game_pooled, Q1_test)>0.5,1,0)

#calculating CER for bad games with pooled model
#Q1_test %>%
#  summarize(CER = mean(Player.Bad.Game != Player.Bad.Pooled.Pred))

#making unpooled model for bad games
Q2_model_bad_game_unpooled <-glmer(Player.Bad.Game~Change.Player.Average.RSI..m.per.s. + (Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_train, family = "binomial")
summary(Q2_model_bad_game_unpooled)

Q1_test$Player.Bad.Unpooled.Pred <- ifelse(predict(Q2_model_bad_game_unpooled, Q1_test, type="response")>0.03,1,0)

Q1_test %>%
  summarize(CER = mean(Player.Bad.Unpooled.Pred != Player.Bad.Game))

Cross Validating Errors for Good and Bad Games with ROC Curves

set.seed(123)

#model to be used for cross validation
model_good_player_pooled_cv <- glm(Player.Good.Game ~ Change.Player.Average.RSI..m.per.s., data=Q1_data, family="binomial")


#Q1_data$Player.Good.Game.Pooled.Pred <- predict(model_good_player_pooled_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Good.Game.Pooled.Pred,
#                   labels = Q1_data$Player.Good.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)

cost <- function(obs, pred){
  mean((pred <=.49) &obs==1 | (pred > 0.49) & obs==0)
}

#cross validating model
good_player_cv <- cv.glm(data=Q1_data,glmfit=model_good_player_pooled_cv,cost,K=5)

#extract average error
good_player_cv$delta[1]

Q1_data$Player.Good.Game.Pooled.Pred <- ifelse(predict(model_good_player_pooled_cv, type="response")>0.49,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Good.Game.Pooled.Pred, y=Q1_data$Player.Good.Game)

set.seed(123)

#model to be used for cross validation
model_good_player_unpooled_cv <- glmer(Player.Good.Game~Change.Player.Average.RSI..m.per.s. + (Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_data, family = "binomial")


#Q1_data$Player.Good.Game.Unpooled.Pred <- predict(model_good_player_unpooled_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Good.Game.Unpooled.Pred,
#                   labels = Q1_data$Player.Good.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)


#making cost function for cross validation for unpooled model
cost <- function(y, yhat){
  mean(((yhat <= 0.46) & y==1) | ((yhat > 0.46) & y==0))
}


#cross validating the unpooled model
summary(cv(model_good_player_unpooled_cv, k=5, clusterVariables = "anon_id", seed=123, criterion=cost))


Q1_data$Player.Good.Game.Unpooled.Pred <- ifelse(predict(model_good_player_unpooled_cv, type="response")>0.46,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Good.Game.Unpooled.Pred, y=Q1_data$Player.Good.Game)

set.seed(123)

#model to be used for cross validation
model_bad_player_unpooled_cv <- glmer(Player.Bad.Game~Change.Player.Average.RSI..m.per.s. + (1|anon_id) + (0+Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_data, family = "binomial")

#Q1_data$Player.Bad.Game.Unpooled.Pred <- predict(model_bad_player_unpooled_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Bad.Game.Unpooled.Pred,
#                   labels = Q1_data$Player.Bad.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)


#making cost function for cross validation for unpooled model
cost <- function(y, yhat){
  mean((yhat <=.05) & y==1 | (yhat > 0.05) & y==0)
}

summary(cv(model_bad_player_unpooled_cv, k=5, clusterVariables = "anon_id", seed=123, criterion=cost))

Q1_data$Player.Bad.Game.Unpooled.Pred <- ifelse(predict(model_bad_player_unpooled_cv, type="response")>0.05,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Bad.Game.Unpooled.Pred, y=Q1_data$Player.Bad.Game)

set.seed(123)

#model to be used for cross validation
model_bad_player_cv <- glm(Player.Bad.Game ~ Change.Player.Average.RSI..m.per.s., data=Q1_data, family="binomial")

#Q1_data$Player.Bad.Game.Pooled.Pred <- predict(model_bad_player_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Bad.Game.Pooled.Pred,
#                   labels = Q1_data$Player.Bad.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)

cost <- function(obs, pred){
  mean((pred <=.03) & obs==1 | (pred > 0.03) & obs==0)
}

#cross validating model
bad_player_cv <- cv.glm(data=Q1_data,glmfit=model_bad_player_cv,cost,K=5)

#extract average error
bad_team_cv$delta[1]

Q1_data$Player.Bad.Game.Pooled.Pred <- if_else(predict(model_bad_player_cv, type="response")>0.03,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Bad.Game.Pooled.Pred, y=Q1_data$Player.Bad.Game)

errors <- data.frame(Model = c("Good-Games-Pooled", "Good-Games-Unpooled", 
                                    "Bad-Games-Pooled", "Bad-Games-Unpooled"),
                     CER = c(0.4503817, 0.3740458   , 0.03053435, 0.2824427),
                     CER.CV = c(0.4837476, 0.456979, 0.2045889, 0.03250478),
                     Classification.Threshold = c(0.49, 0.46, 0.03, 0.05))

kable(errors)

Based on the errors that come from the cross validation of all of the models, we can see that the unpooled models outperform the pooled models. By fitting a unique model to each player we were able to get better predictions for in game success than fitting a common model to all of the players. This suggests that RSI may have differing predictive power for in game performance for each player. This means that changes in RSI may be a more powerful predictor for in game success than for other athletes on the team. This also suggests that some information may be lost when averaging changes in RSI for the team. One thing to note about these models is that the metrics shown above are predicted rates in which the model will make a wrong prediction. It’s important to know that each model gives the predicted probability that the given observation is a “good game” for the first two models or the probability that a given observation is a “bad game” for the last two models. Normally the threshold to consider a prediction to be a 1 or a 0 is if the probability is above or below 0.5 (rounded to the nearest whole number). These models were given different classification thresholds. This was for computational reasons since “bad games” for players are very rare and will therefore be harder to identify. For the models that are predicting whether a game will be a “good game” or not, they don’t have a lot of predictive power.

Question 3:

Is the previous week’s load related to RSI?

Recalculating Gameday Minus Columns

# Recalculating the Gameday.Minus.N.Acceleration.Load Columns

# Create a df with relevant columns and add column specifying if that day was a game or not
daily_load <- kinexon_session_clean |>
  mutate(Game.Day = if_else(Type == "Game", 1, 0)) |>
  filter(!is.na(Daily.Total.Accel.Load.Accum)) |>
  select(anon_id, Date, Daily.Total.Accel.Load.Accum, Game.Day) |>
  arrange(anon_id, Date)

# Get all game day dates per player as a lookup table
game_days <- daily_load |>
  filter(Game.Day == 1) |>
  select(anon_id, Game.Date = Date)
    

# Add column Day.Type for GD-1 - GD-4 and gamedays
daily_load <- daily_load |>
  left_join(game_days, by = "anon_id", relationship = "many-to-many") |>
  mutate(
    Days.Until.Game = as.numeric(difftime(Game.Date, Date, units = "days")),
    Day.Type = case_when(
      Days.Until.Game == 1 ~ "GD-1",
      Days.Until.Game == 2 ~ "GD-2",
      Days.Until.Game == 3 ~ "GD-3",
      Days.Until.Game == 4 ~ "GD-4",
      Game.Day == 1 ~ "GameDay",
      TRUE ~ NA_character_
    )
  ) |>
  filter(!is.na(Day.Type)) |>    # keep only rows within 4 days before game
  select(-Game.Date, -Days.Until.Game)

# Add columns for Gameday.Minus.N.Acceleration.Load
Gameday_Minus <- daily_load |>
  mutate(
    Gameday.Minus.1.Acceleration.Load = if_else(Day.Type == "GD-1", Daily.Total.Accel.Load.Accum, NA_real_),
    Gameday.Minus.2.Acceleration.Load = if_else(Day.Type == "GD-2", Daily.Total.Accel.Load.Accum, NA_real_),
    Gameday.Minus.3.Acceleration.Load = if_else(Day.Type == "GD-3", Daily.Total.Accel.Load.Accum, NA_real_),
    Gameday.Minus.4.Acceleration.Load = if_else(Day.Type == "GD-4", Daily.Total.Accel.Load.Accum, NA_real_),
    GameDay.Acceleration.Load = if_else(Day.Type == "GameDay", Daily.Total.Accel.Load.Accum, NA_real_)
  )

# Delete duplicate rows
Gameday_Minus <- Gameday_Minus |>
  distinct()

# Create columns for Average.Gameday.Minus.N.Acceleration.Load

# Calculate and add columns for running averages
Gameday_Minus <- Gameday_Minus |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  mutate(
    # Cumulative count and sum for GD-1
    cum_count_gd1 = cumsum(!is.na(Gameday.Minus.1.Acceleration.Load)),
    cum_sum_gd1 = cumsum(replace_na(Gameday.Minus.1.Acceleration.Load, 0)),
    Average.Gameday.Minus.1.Acceleration.Load = if_else(
      cum_count_gd1 > 0,
      cum_sum_gd1 / cum_count_gd1,
      NA_real_
    ),
    
    # Cumulative count and sum for GD-2
    cum_count_gd2 = cumsum(!is.na(Gameday.Minus.2.Acceleration.Load)),
    cum_sum_gd2 = cumsum(replace_na(Gameday.Minus.2.Acceleration.Load, 0)),
    Average.Gameday.Minus.2.Acceleration.Load = if_else(
      cum_count_gd2 > 0,
      cum_sum_gd2 / cum_count_gd2,
      NA_real_
    ),
    
    # Cumulative count and sum for GD-3
    cum_count_gd3 = cumsum(!is.na(Gameday.Minus.3.Acceleration.Load)),
    cum_sum_gd3 = cumsum(replace_na(Gameday.Minus.3.Acceleration.Load, 0)),
    Average.Gameday.Minus.3.Acceleration.Load = if_else(
      cum_count_gd3 > 0,
      cum_sum_gd3 / cum_count_gd3,
      NA_real_
    ),
    
    # Cumulative count and sum for GD-4
    cum_count_gd4 = cumsum(!is.na(Gameday.Minus.4.Acceleration.Load)),
    cum_sum_gd4 = cumsum(replace_na(Gameday.Minus.4.Acceleration.Load, 0)),
    Average.Gameday.Minus.4.Acceleration.Load = if_else(
      cum_count_gd4 > 0,
      cum_sum_gd4 / cum_count_gd4,
      NA_real_
    )
  ) |>
  ungroup() |>
  select(-starts_with("cum_count_"), -starts_with("cum_sum_")) # clean up columns

# Calculate running averages of Gameday.Acceleration.Load and add column Average.Gameday.Acceleration.Load

Gameday_Minus <- Gameday_Minus |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  mutate(
    cum_count_gd = cumsum(!is.na(GameDay.Acceleration.Load)),
    cum_sum_gd = cumsum(replace_na(GameDay.Acceleration.Load, 0)),
    Average.GameDay.Acceleration.Load = if_else(
      cum_count_gd > 0,
      cum_sum_gd / cum_count_gd,
      NA_real_
    )
  ) |>
  ungroup() |>
  select(-cum_count_gd, -cum_sum_gd) # clean up columns

Calculate load of 7, 3, and 1, day before test day including test date.

# Filter to only have data from this season
RSI <- VALD_ForceDecks_clean |>
  filter(Date >= as.Date("2024-10-10")) |>
  select(anon_id, Date, RSI..m.per.s., Trial) |>
  arrange(anon_id, Date) |>
  group_by(anon_id, Date) |>
  mutate(Daily.Avg.RSI..m.per.s = mean(RSI..m.per.s.)) |>
  ungroup()


# Plot RSI over season
ggplot(RSI, aes(x=Date, y=Daily.Avg.RSI..m.per.s, color=anon_id)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  labs(title = "Athlete RSI Scores Throughout Season") +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# Create load data frame containing columns relevant to load
load <- kinexon_session_clean |>
  filter(Date >= as.Date("2024-10-10")) |>
  select(anon_id, Date, Daily.Total.Accel.Load.Accum, Jump.Load) |>
  filter(!is.na(Daily.Total.Accel.Load.Accum))


# Look into Daily Acceleration Load for one athlete
id_40_load <- load |>
  filter(anon_id == "ID_40")
# Plot of ID_40's Daily Accel Load
ggplot(data = id_40_load, aes(x = Date, y = Daily.Total.Accel.Load.Accum)) +
  geom_point() +
  geom_line() +
  labs(title = "ID_40 Daily Acceleration Load") +
  theme_classic()

# Plot of all athlete's daily accel load throughout season
ggplot(data = load, aes(x=Date, y=Daily.Total.Accel.Load.Accum, color=anon_id)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "2024-25 Baskteball Athlete's Daily Acceleration Load") +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# Plot of all athlete's daily jump load throughout season
ggplot(data = load, aes(x=Date, y=Jump.Load, color=anon_id)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "2024-25 Baskteball Athlete's Daily Jump Load") +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

The graphs seem to show that for most athletes, RSI scores increased throughout the season, while daily acceleration load decreased and jump load about stayed the same.

# Calculate the average and total acceleration and jump load in the last 7, 3, and 2 days.

## Last 7 days not last 7 sessions
# Include that day in the last 7. (That day and the 6 days before that day, total of 7)
rolling_NDay_load <- load |>
  filter(!is.na(Daily.Total.Accel.Load.Accum)) |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  group_split() |>
  map_dfr(~{
    athlete_data <- .
    map_dfr(seq_len(nrow(athlete_data)), function(i) {
      current_row <- athlete_data[i, ]
      past_week <- athlete_data |>
        filter(Date <= current_row$Date & Date > current_row$Date - days(7)) # (<=) includes that day in the last 7.
      
      current_row |>
        mutate(
          Avg..7Day..Accel.Load = mean(past_week$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Total..7Day..Accel.Load = sum(past_week$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Avg..7Day..Jump.Load = mean(past_week$Jump.Load, na.rm = TRUE),
          Total..7Day..Jump.Load = sum(past_week$Jump.Load, na.rm = TRUE),
          Days.Used.7Day = nrow(past_week)
        )
    })
  })



# Add 3 Day Load
# That day and the two days before it
rolling_NDay_load <- rolling_NDay_load |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  group_split() |>
  map_dfr(~{
    athlete_data <- .
    map_dfr(seq_len(nrow(athlete_data)), function(i) {
      current_row <- athlete_data[i, ]
      past_3days <- athlete_data |>
        filter(Date <= current_row$Date & Date > current_row$Date - days(3))  # includes current day

      current_row |>
        mutate(
          Avg..3Day..Accel.Load = mean(past_3days$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Total..3Day..Accel.Load = sum(past_3days$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Avg..3Day..Jump.Load = mean(past_3days$Jump.Load, na.rm = TRUE),
          Total..3Day..Jump.Load = sum(past_3days$Jump.Load, na.rm = TRUE),
          Days.Used.3Day = nrow(past_3days)
        )
    })
  })



# Add 2 Day Load
# That day and the day before it
rolling_NDay_load <- rolling_NDay_load |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  group_split() |>
  map_dfr(~{
    athlete_data <- .
    map_dfr(seq_len(nrow(athlete_data)), function(i) {
      current_row <- athlete_data[i, ]
      past_2days <- athlete_data |>
        filter(Date <= current_row$Date & Date > current_row$Date - days(2))  # includes current day

      current_row |>
        mutate(
          Avg..2Day..Accel.Load = mean(past_2days$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Total..2Day..Accel.Load = sum(past_2days$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Avg..2Day..Jump.Load = mean(past_2days$Jump.Load, na.rm = TRUE),
          Total..2Day..Jump.Load = sum(past_2days$Jump.Load, na.rm = TRUE),
          Days.Used.2Day = nrow(past_2days)
        )
    })
  })

RSI_daily_avg <- RSI |>
  distinct(anon_id, Date, Daily.Avg.RSI..m.per.s)

# Merge the RSI and rolling day load data frames
load_and_RSI <- RSI_daily_avg |>
  left_join(
    rolling_NDay_load,
    by = "anon_id",
    relationship = "many-to-many"
  ) |>
  filter(Date.y <= Date.x & Date.y >= Date.x - days(7)) |> # <= includes same day. On days when RSI and accel load are measured, accel load is measured first, meaning we should include it.
  group_by(anon_id, Date.x) |>
  slice_max(Date.y, n = 1) |>  # Keep the most recent load data before RSI
  ungroup() |>
  rename(
    RSI_Date = Date.x,
    Load_Date = Date.y,
    RSI = Daily.Avg.RSI..m.per.s
  )


# Now that only have load data for N days before a test (and that day) we rename the columns
load_and_RSI <- load_and_RSI |> rename(
  `TD-6..Avg..Accel.Load` = Avg..7Day..Accel.Load,
  `TD-6..Total..Accel.Load` = Total..7Day..Accel.Load,
  `TD-6..Avg..Jump.Load` = Avg..7Day..Jump.Load,
  `TD-6..Total..Jump.Load` = Total..7Day..Jump.Load,

  `TD-2..Avg..Accel.Load` = Avg..3Day..Accel.Load,
  `TD-2..Total..Accel.Load` = Total..3Day..Accel.Load,
  `TD-2..Avg..Jump.Load` = Avg..3Day..Jump.Load,
  `TD-2..Total..Jump.Load` = Total..3Day..Jump.Load,

  `TD-1..Avg..Accel.Load` = Avg..2Day..Accel.Load,
  `TD-1..Total..Accel.Load` = Total..2Day..Accel.Load,
  `TD-1..Avg..Jump.Load` = Avg..2Day..Jump.Load,
  `TD-1..Total..Jump.Load` = Total..2Day..Jump.Load
)

# TD-6 means the load data from the last 6 days before the test day and the test day
# (load data from the last 6 days and the test day)
# TD-2 means the load data from the two days prior to the test day and the test day
# TD-1 means the data from the day before the test day and data from the test day

# Correlation Tests for each variable vs RSI

# Define target variables
load_vars <- c(
  "TD-6..Avg..Accel.Load", "TD-6..Total..Accel.Load",
  "TD-2..Avg..Accel.Load", "TD-2..Total..Accel.Load",
  "TD-1..Avg..Accel.Load", "TD-1..Total..Accel.Load",
  "TD-6..Avg..Jump.Load",  "TD-6..Total..Jump.Load",
  "TD-2..Avg..Jump.Load",  "TD-2..Total..Jump.Load",
  "TD-1..Avg..Jump.Load",  "TD-1..Total..Jump.Load"
)

# Compute correlations
correlation_results <- map_dfr(load_vars, function(var) {
  pearson <- cor(load_and_RSI$RSI, load_and_RSI[[var]], use = "complete.obs")
  spearman <- cor(load_and_RSI$RSI, load_and_RSI[[var]], method = "spearman", use = "complete.obs")
  
  tibble(
    Variable = var,
    Pearson = round(pearson, 3),
    Spearman = round(spearman, 3)
  )
})


kable(correlation_results)

Variable	Pearson	Spearman
TD-6..Avg..Accel.Load	0.008	0.023
TD-6..Total..Accel.Load	0.042	0.055
TD-2..Avg..Accel.Load	0.019	0.027
TD-2..Total..Accel.Load	-0.014	-0.042
TD-1..Avg..Accel.Load	-0.048	-0.047
TD-1..Total..Accel.Load	-0.025	-0.010
TD-6..Avg..Jump.Load	0.258	0.231
TD-6..Total..Jump.Load	0.238	0.230
TD-2..Avg..Jump.Load	0.243	0.229
TD-2..Total..Jump.Load	0.160	0.105
TD-1..Avg..Jump.Load	0.210	0.205
TD-1..Total..Jump.Load	0.189	0.180

NA

# AVERAGE ACCEL LOAD:
# Relationship between RSI and 7-Day Average Acceleration Load
ggplot(load_and_RSI, aes(x = `TD-6..Avg..Accel.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 7-Day Avg Accel Load",
    x = "7-Day Avg Accel Load (TD-6 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# Relationship between RSI and 3-Day Average Acceleration Load
ggplot(load_and_RSI, aes(x = `TD-2..Avg..Accel.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 3-Day Avg Accel Load",
    x = "3-Day Avg Accel Load (TD-2 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# Relationship between RSI and 2-Day Average Acceleration Load
ggplot(load_and_RSI, aes(x = `TD-1..Avg..Accel.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 2-Day Avg Accel Load",
    x = "2-Day Avg Accel Load (TD-1 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# TOTAL ACCEL LOAD:
# RSI and 7-Day Total Acceleration Load
ggplot(load_and_RSI, aes(x = `TD-6..Total..Accel.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 7-Day Avg Accel Load",
    x = "7-Day Total Accel Load (TD-6 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# RSI and 3-Day Total Acceleration Load
ggplot(load_and_RSI, aes(x = `TD-2..Total..Accel.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 3-Day Avg Accel Load",
    x = "3-Day Total Accel Load (TD-2 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# RSI and 2-Day Average Acceleration Load
ggplot(load_and_RSI, aes(x = `TD-1..Total..Accel.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 2-Day Avg Accel Load",
    x = "2-Day Total Accel Load (TD-1 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# AVERAGE JUMP LOAD:
# RSI and 7-Day Average Jump Load
ggplot(load_and_RSI, aes(x = `TD-6..Avg..Jump.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 7-Day Avg Jump Load",
    x = "7-Day Avg Jump Load (TD-6 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# RSI and 3-Day Average Jump Load
ggplot(load_and_RSI, aes(x = `TD-2..Avg..Jump.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 3-Day Avg Jump Load",
    x = "3-Day Avg Jump Load (TD-2 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# RSI and 2-Day Average Jump Load
ggplot(load_and_RSI, aes(x = `TD-1..Avg..Jump.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 2-Day Avg Jump Load",
    x = "2-Day Avg Jump Load (TD-1 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# TOTAL JUMP LOAD:
# RSI and 7-Day Total Jump Load
ggplot(load_and_RSI, aes(x = `TD-6..Total..Jump.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 7-Day Avg Jump Load",
    x = "7-Day Total Jump Load (TD-6 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# RSI and 3-Day Total Jump Load
ggplot(load_and_RSI, aes(x = `TD-2..Total..Jump.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 3-Day Avg Jump Load",
    x = "3-Day Total Jump Load (TD-2 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# RSI and 2-Day Total Jump Load
ggplot(load_and_RSI, aes(x = `TD-1..Total..Jump.Load`, y = RSI)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE, color = "#CFB87C") +
  labs(
    title = "RSI vs. 2-Day Avg Jump Load",
    x = "2-Day Total Jump Load (TD-1 to TD)",
    y = "RSI (m/s)"
  ) +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# Plot for each athlete
unique_ids <- unique(load_and_RSI$anon_id)

for (id in unique_ids) {
  p <- ggplot(filter(load_and_RSI, anon_id == id),
              aes(x = `TD-6..Avg..Accel.Load`, y = RSI)) +
    geom_point() +
    geom_smooth(method = "loess", color = "#CFB87C") +
    labs(title = paste("RSI vs. Prior Week's Accel Load:", id),
         x = "Total Load Last 7 Days",
         y = "RSI")
   print(p)
   
}

`geom_smooth()` using formula = 'y ~ x'

Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  span too small.   fewer data values than degrees of freedom.
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at 710.48
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 133.28
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  0
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 6076.6
Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata)) as.matrix(model.frame(delete.response(terms(object)),  :
  span too small.   fewer data values than degrees of freedom.
Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata)) as.matrix(model.frame(delete.response(terms(object)),  :
  pseudoinverse used at 710.48
Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata)) as.matrix(model.frame(delete.response(terms(object)),  :
  neighborhood radius 133.28
Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata)) as.matrix(model.frame(delete.response(terms(object)),  :
  reciprocal condition number  0
Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x else if (is.data.frame(newdata)) as.matrix(model.frame(delete.response(terms(object)),  :
  There are other near singularities as well. 6076.6
Warning in max(ids, na.rm = TRUE) :
  no non-missing arguments to max; returning -Inf

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  Chernobyl! trL>n 6
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  Chernobyl! trL>n 6
Warning in sqrt(sum.squares/one.delta) : NaNs produced
Warning in stats::qt(level/2 + 0.5, pred$df) : NaNs produced
Warning in max(ids, na.rm = TRUE) :
  no non-missing arguments to max; returning -Inf

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

# Correlation Tests for each variable vs RSI

# Calculate overall Pearson correlations for each variable
overall_corrs <- map_dfr(load_vars, ~ {
  load_and_RSI |> 
    summarize(
      r = cor(RSI, .data[[.x]], use = "complete.obs"),
      n = n()
    ) |> 
    mutate(variable = .x)
})

# Basic plot (no theming)
ggplot(overall_corrs, aes(x = reorder(variable, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
            hjust = ifelse(overall_corrs$r >= 0, -0.1, 1.1),
            size = 3) +
  coord_flip() +
  labs(title = "Correlation Between Load Metrics and RSI",
       x = "Load Metric",
       y = "Pearson Correlation (r)") +
  scale_x_discrete(labels = function(x) gsub("\\.\\.", ": ", x)) +  # Clean names
  theme_classic()

NA
NA

Jump Load seems to be more correlated and even positively correlated than accel load.

# Per athlete accel load correlations with RSI

# ACCEL LOAD:
# Average 7 day accel load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-6..Avg..Accel.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_avg$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Avg 7-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") + 
  theme_classic()

# Total 7 day accel load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-6..Total..Accel.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_total$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Total 7-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Average 3 day accel load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-2..Avg..Accel.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_avg$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Avg 3-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Total 3 day accel load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-2..Total..Accel.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_total$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Total 3-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Average 2 day accel load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-1..Avg..Accel.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_avg$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Avg 2-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Total 2 day accel load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-1..Total..Accel.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_total$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Total 2-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

ID_64 seem to have a moderatley strong negative correlation between acceleration load and RSI.

# JUMP LOAD
# Average 7 day jump load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-6..Avg..Jump.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_avg$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Avg 7-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Total 7 day jump load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-6..Total..Jump.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_total$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Total 7-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Average 3 day jump load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-2..Avg..Jump.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_avg$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Avg 3-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Total 3 day jump load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-2..Total..Jump.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_total$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Total 3-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Average 2 day jump load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-1..Avg..Jump.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_avg$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Avg 2-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# Total 2 day jump load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, `TD-1..Total..Jump.Load`, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "#CFB87C") +
  geom_text(aes(label = round(r, 2)), 
              hjust = ifelse(per_athlete_corrs_total$r >= 0, -0.1, 1.1),
              size = 3.2,
              color = "#000000") +
  coord_flip() +
  labs(title = "Per Athlete Correlation Between Total 2-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)") +
  theme_classic()

# This allows us to see which athlete's RSI scores are most effected by the prior week's load

set.seed(42)

# Create grouped 5-fold CV splits by anon_id
cv_splits <- group_vfold_cv(load_and_RSI, group = anon_id, v = 5)

# Mixed Effects Model with standardized predictors
cv_mixed <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  model <- lmer(RSI ~ z_avg_accel + z_avg_jump + (1 | anon_id), data = train)
  preds <- predict(model, newdata = test, allow.new.levels = TRUE)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# GAM Model with standardized predictors
cv_gam <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  model <- gam(RSI ~ s(z_avg_accel) + s(z_avg_jump), data = train)
  preds <- predict(model, newdata = test)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# Run CV for both models
mixed_rmses <- map_dbl(cv_splits$splits, cv_mixed)

Warning: Unknown or uninitialised column: `Last.Week..Avg..Accel.Load`.
Warning in mean.default(train$Last.Week..Avg..Accel.Load, na.rm = TRUE) :
  argument is not numeric or logical: returning NA
Warning: Unknown or uninitialised column: `Last.Week..Avg..Accel.Load`.
Warning: Unknown or uninitialised column: `Last.Week..Avg..Jump.Load`.
Warning in mean.default(train$Last.Week..Avg..Jump.Load, na.rm = TRUE) :
  argument is not numeric or logical: returning NA
Warning: Unknown or uninitialised column: `Last.Week..Avg..Jump.Load`.

Error in `map_dbl()`:
ℹ In index: 1.
Caused by error in `mutate()`:
ℹ In argument: `z_avg_accel = (Last.Week..Avg..Accel.Load -
  mean_accel)/sd_accel`.
Caused by error:
! object 'Last.Week..Avg..Accel.Load' not found
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.

# Fit fixed effects model on full data
avg_fixed_model <- lm(RSI ~ Last.Week..Avg..Accel.Load + Last.Week..Avg..Jump.Load + factor(anon_id),
                  data = load_and_RSI)

Error in eval(predvars, data, env) : 
  object 'Last.Week..Avg..Accel.Load' not found

set.seed(42)

# Create grouped 5-fold CV splits by anon_id
cv_splits <- group_vfold_cv(load_and_RSI, group = anon_id, v = 5)

# Mixed Effects Model with standardized total load predictors
cv_mixed_total <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Total..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Total..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Total..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Total..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  model <- lmer(RSI ~ z_total_accel + z_total_jump + (1 | anon_id), data = train)
  preds <- predict(model, newdata = test, allow.new.levels = TRUE)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# GAM Model with standardized total load predictors
cv_gam_total <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Total..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Total..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Total..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Total..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  model <- gam(RSI ~ s(z_total_accel) + s(z_total_jump), data = train)
  preds <- predict(model, newdata = test)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# Run CV for both models with total loads
mixed_rmses_total <- map_dbl(cv_splits$splits, cv_mixed_total)

Warning: Unknown or uninitialised column: `Last.Week..Total..Accel.Load`.
Warning in mean.default(train$Last.Week..Total..Accel.Load, na.rm = TRUE) :
  argument is not numeric or logical: returning NA
Warning: Unknown or uninitialised column: `Last.Week..Total..Accel.Load`.
Warning: Unknown or uninitialised column: `Last.Week..Total..Jump.Load`.
Warning in mean.default(train$Last.Week..Total..Jump.Load, na.rm = TRUE) :
  argument is not numeric or logical: returning NA
Warning: Unknown or uninitialised column: `Last.Week..Total..Jump.Load`.

Error in `map_dbl()`:
ℹ In index: 1.
Caused by error in `mutate()`:
ℹ In argument: `z_total_accel = (Last.Week..Total..Accel.Load -
  mean_accel)/sd_accel`.
Caused by error:
! object 'Last.Week..Total..Accel.Load' not found
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.

Question 4

What is each athlete’s variation in RSI? What is a meaningful change in RSI for the team and for individual athletes?

# Some basic visualizations

VALD_ForceDecks_clean <- VALD_ForceDecks_clean %>%
  filter(Date >= as.Date("2024-09-20", '%Y-%m-%d')) %>%
  filter(!(anon_id %in% c("ID_10", "ID_14", "ID_15", "ID_18", "ID_54", "ID_64")))
  
# Boxplot showing RSI distribution for each athlete
ggplot(VALD_ForceDecks_clean, aes(anon_id, RSI..m.per.s.)) +
  geom_boxplot(color = "#CFB87C") +
  labs(title = "RSI Variability by Athlete") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme_classic()

# Line and scatter plot showing RSI trends over time for each athlete
# Includes a black line for the team average RSI to compare individual performance to
ggplot(VALD_ForceDecks_clean, aes(Date, Player.Average.RSI..m.per.s., color = anon_id)) +
  geom_point() + 
  geom_smooth(method = "lm", se=FALSE) +
  labs(title = "Player RSI Over Time with Team Average Line") +
  geom_line(aes(y = Team.Average.RSI..m.per.s.), 
            color = "black") +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

# Smoothed RSI Trend over time
ggplot(load_and_RSI, aes(x = RSI_Date, y = RSI)) +
  geom_smooth(method = "loess", color = "#CFB87C", fill = "#A2A4A3") +
  labs(title = "Smoothed RSI Trend",
       x = "Date",
       y = "RSI Value") +
  theme_classic()

`geom_smooth()` using formula = 'y ~ x'

kendall_cor = cor(VALD_ForceDecks_clean$X, VALD_ForceDecks_clean$Team.Average.RSI..m.per.s., method = "kendall")*-1
ggplot(data=VALD_ForceDecks_clean, aes(x=Date, y=Player.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_smooth(method="lm") +
  labs(title=kendall_cor)

`geom_smooth()` using formula = 'y ~ x'

#taking out all of the IDs in the data set
IDs <- unique(VALD_ForceDecks_clean$anon_id)

#for all unique IDs, calculate Kendall rank correlation, plot player RSI throughout season 
for(i in 1:11){
  RSI <- unique(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.)
  kendall_cor <- cor(1:length(RSI), RSI, method="kendall")*-1 #calculating Kendall correlation

#plotting each player's RSI throughout the season with their correlation and rough model of the trend
  p <- ggplot(data = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],], aes(x=Date, y=Player.Average.RSI..m.per.s.)) + 
    geom_point() +
    geom_smooth(method="lm") +
    labs(title = kendall_cor, subtitle = IDs[i])
  
#pring out each plot
  print(p)
}

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

#for all the players, calculate mean and variance of their RSI measurements and plot
for(i in 1:16){
  #for every player calculate mean and variance using time series assumptions and functions from the time series functions file
  mv <- sample_mean_and_variance(unique(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.))
  
  #plotting RSI measurements along with calculated mean and 1 standard deviation above and below mean for each player along with their measured RSIs throughout the season
  p <- ggplot(data = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],], aes(x=Date, y=Player.Average.RSI..m.per.s.)) + 
    geom_point() +
    geom_hline(yintercept = mv[[1]]) +
    geom_hline(yintercept = mv[[1]] + sqrt(mv[[2]]), color="#CFB87C") +
    geom_hline(yintercept = mv[[1]] - sqrt(mv[[2]]), color="#CFB87C") +
    labs(title = "RSI Measurements Throughout Season with Mean and 1 Standard Deviation",subtitle = IDs[i])
  
  #printing out each plot
  print(p)
}

#making empty columns to be filled with indicator variables
VALD_ForceDecks_clean$Player.Above.1.SD <- rep(NA,565)
VALD_ForceDecks_clean$Player.Below.1.SD <- rep(NA,565)
VALD_ForceDecks_clean$Player.Above.2.SD <- rep(NA,565)
VALD_ForceDecks_clean$Player.Below.2.SD <- rep(NA,565)

#for every athlete, checking if a specific RSI measurement is outside 1 standard deviation or 2 calculated for each athlete individually. 
for(i in 1:16){
  #calculating each player's mean and variance with time series function on the time series functions file
  mv <- sample_mean_and_variance(unique(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.))
  
  #calculating thresholds for 1 and 2 standard deviations away from mean for each player
  sdu1 <- mv[[1]] + sqrt(mv[[2]])
  sdl1 <- mv[[1]] - sqrt(mv[[2]])
  sdu2 <- mv[[1]] + (2*sqrt(mv[[2]]))
  sdl2 <- mv[[1]] - (2*sqrt(mv[[2]]))
  
#checking if RSI value is below 1 standard deviation for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.1.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.<sdl1), 1, 0)
  
#checking if RSI value is above 1 standard deviation for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.1.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.>sdu1), 1, 0)
  
#checking if RSI value is below 2 standard deviations for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.2.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.<sdl2), 1, 0)
  
#checking if RSI value is above 2 standard deviations for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.2.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.>sdu2), 1, 0)
}

#making a plot to show potentially significant RSI measurements throughout the season
#dark green means point is above 1 standard deviation from mean
#light green means point is above 2 standard deviations from mean
#dark red means point is below 1 standard deviation from mean
#light red means point is below 2 standard deviations from mean
ggplot(data=VALD_ForceDecks_clean, aes(x=Date, Player.Average.RSI..m.per.s.)) +
  geom_point(color = ifelse(VALD_ForceDecks_clean$Player.Above.2.SD == 1, "palegreen", ifelse(VALD_ForceDecks_clean$Player.Below.2.SD == 1, "red", ifelse(VALD_ForceDecks_clean$Player.Below.1.SD == 1, "red4", ifelse(VALD_ForceDecks_clean$Player.Above.1.SD == 1, "palegreen4", "grey80")))), alpha = 0.4) +
  labs(title = "Significant RSI Measurements Throughout the Season") +
  theme_minimal()

This plot shows each player’s measured RSI throughout the season. Points that are colored gray are points where the player did not have a significant difference in their RSI from their mean RSI throughout the season. Dark green dots show RSI measurements that are above 1 standard deviation for each player and light green dots are RSI measurements that are above 2 standard deviations for each player. Along the same lines, dark red dots shown points in which the player had an RSI below 1 standard deviation from their RSI distribution and light red dots show points in which the player’s measured RSI is below 2 standard deviations of their RSI distribution. One interesting thing about this plot is that there are no measurements that are below even 1 standard deviation for any player after February starts. If anything, there seems to be a higher concentration of RSI measurements that are at least above 1 standard deviation from the mean in the last 5 measurements of the season, after February starts than in the rest of the season. This plot though, does not help us understand how players relate to each other. This chart along with the box plots of each player’s RSI measurements support the conclusion that each player has a very unique RSI distribution. The plot above gets broken down by each player below.

#cycling through all of the players and plotting their RSI measurements throughout the season
for(i in 1:11){
  #plot the measured RSI for each player throughout the season
  #dark green means point is above 1 standard deviation from mean
#light green means point is above 2 standard deviations from mean
#dark red means point is below 1 standard deviation from mean
#light red means point is below 2 standard deviations from mean
p <- ggplot(data=VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],], aes(x=Date, Player.Average.RSI..m.per.s.)) +
  geom_point(color = ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.2.SD == 1, "palegreen", ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.2.SD == 1, "red", ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.1.SD == 1, "red4", ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.1.SD == 1, "palegreen4", "grey70"))))) +
  ylim(0.7311404, 2.963925) +
  labs(title = "Measured RSI Values Throughout the Season and Their Significance", subtitle = IDs[i]) +
  theme_minimal()

print(p)
}

The plots above help us to see that the RSI distributions for each player are going to be wildly different from each other. In fact, some of the plots have no overlap with each other (e.g. ID_5 and ID_74). This suggests that defining one singular threshold for what is a meaningful change in RSI will not be accurate for all of the players. These plots when looked at all together also help us to see that when players have an RSI measurement that is below one standard deviation from their mean, it is usually associated with another measurement that is below one standard deviation either directly or in the following two or three measurements before returning to a measurement that is close to their season mean. This is not true for all players but seems to be true for most of the players. This suggests that for most of the team the recovery time for neuromuscular fatigue that results in a noticeably lower RSI score, might be around 8 to 10 days given that RSI measurements are taken every 4 to 5 days on average.

VALD_ForceDecks_clean <- VALD_ForceDecks_clean %>%
  group_by(Date) %>%
  #making 4 new columns that will count how many players were in each category for each date
  mutate(Num.Above.1SD = round(mean(Player.Above.1.SD)*length(unique(anon_id))),
         Num.Below.1SD = round(mean(Player.Below.1.SD)*length(unique(anon_id))),
         Num.Above.2SD = round(mean(Player.Above.2.SD)*length(unique(anon_id))),
         Num.Below.2SD = round(mean(Player.Below.2.SD)*length(unique(anon_id)))) %>%
  ungroup()

Dates <- unique(VALD_ForceDecks_clean$Date)
rows <- rep(NA, 22)

for(i in 1:22){
  rows[i] <- min(VALD_ForceDecks_clean[VALD_ForceDecks_clean$Date == Dates[i],]$X)
}

test <- VALD_ForceDecks_clean[VALD_ForceDecks_clean$X %in% rows,]

ggplot(data=test, aes(x=Date)) +
  geom_col(aes(y=Num.Above.1SD), fill = "palegreen4") +
  geom_col(aes(y=Num.Above.2SD), fill = "palegreen") +
  ylim(0,7) +
  labs(title = "Number of Players Above Their Mean Season RSI Measurement",
       y="Number Above Mean RSI Measurement") +
  theme_minimal()

ggplot(data=test, aes(x=Date)) +
  geom_col(aes(y=Num.Below.1SD), fill = "red4") +
  geom_col(aes(y=Num.Below.2SD), fill = "red") +
  ylim(0,7)+
  labs(title = "Number of Players Below Their Mean Season RSI Measurement",
       y="Number Below Mean RSI Measurement") +
  theme_minimal()

From the plots above, we can see that there are two general parts of the season in which the team has more frequent drops in RSI and the team has more frequent gains in RSI. These plots make it easy to see that the team had a lot more frequent drops in RSI in the beginning of the season than at the end of the season. In fact, the last major drop for any player on the team happened at the end of January and there are no more drops for the rest of the season. With this, we see that RSI gains are not as common throughout the beginning of the season. We only see one or two players that have an RSI measurement above one or two standard deviations above their mean until February. It is only once February starts that we see an increase in the number of players that have an RSI at least one standard deviation above their mean.

Looking at player average RSI values

#calculating the mean and variance of RSI for the players throughout the season
team_mean <- mean(VALD_ForceDecks_clean$Player.Average.RSI..m.per.s.)
team_var <- sample_mean_and_variance(VALD_ForceDecks_clean$Player.Average.RSI..m.per.s.)[[2]]

#plotting RSI measurements throughout the season, adding in the mean RSI value and 1 standard deviation above and below the mean
ggplot(VALD_ForceDecks_clean, aes(x=Date, y=Player.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = team_mean) +
  geom_hline(yintercept = team_mean - sqrt(team_var), color="#CFB87C") +
  geom_hline(yintercept = team_mean + sqrt(team_var), color="#CFB87C")

#plotting the RSI values throughout the season and plotting the team RSI trend overtop to see if there are underying trends that could be explained with the team trend
ggplot(VALD_ForceDecks_clean, aes(x=Date, y=Player.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_line(data=VALD_ForceDecks_clean, aes(Date,Team.Average.RSI..m.per.s.), color = "#CFB87C")

#cycle through all of the players
for(i in 1:11){
  #making a plot with date on x axis and average 7 day jump load on the y axis to see if there is potentially any relationship in the trend of the load and changes in RSI
r <- ggplot(rolling_NDay_load[rolling_NDay_load$anon_id == IDs[i],], aes(x=Date, y=Total..7Day..Jump.Load)) +
  geom_point() +
  ylim(0, 400000) +
  labs(title="Jump Loads through the Season", subtitle = IDs[i]) +
  #calculating a moving average to better see the underlying shape of the data
  geom_line(y=moving_average(rolling_NDay_load[rolling_NDay_load$anon_id == IDs[i],]$Total..7Day..Jump.Load,3),
            color = "#CFB87C") +
  
  #adding in vertical lines that represent the dates in which the player had an RSI measurement above or below their mean by at least 1 SD. Same color palette applies
  geom_vline(xintercept = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id==IDs[i] & VALD_ForceDecks_clean$Player.Above.1.SD==1,]$Date, color = "palegreen4") +
  geom_vline(xintercept = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id==IDs[i] & VALD_ForceDecks_clean$Player.Above.2.SD==1,]$Date, color = "palegreen") + 
  geom_vline(xintercept = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id==IDs[i] & VALD_ForceDecks_clean$Player.Below.1.SD==1,]$Date, color = "red4") +
  geom_vline(xintercept = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id==IDs[i] & VALD_ForceDecks_clean$Player.Below.2.SD==1,]$Date, color = "red") +
  geom_hline(yintercept = mean(rolling_NDay_load[rolling_NDay_load$anon_id == IDs[i],]$Total..7Day..Jump.Load))

print(r)
}

Looking at the plots above, there does not seem to be a clear relationship between RSI measurements and loading summarized by any length (2,3,7 days). This suggests that there may have been some other change in training that had a positive impact on RSI measurements for the entire team. This may be due to something that is not loading related. Blake said it might be more targeted lifts.